In [54]:
# load in template time series to replace data
templateTS=nb.load('C:/Users/adam/template-rest_p2mm_masked.dtseries.nii')
templateData=templateTS.get_fdata()
# extract cortex
templateDataCort=templateData[:,hcp.struct.cortex]
# sampling frequency
TR=.8
fs=(1/TR)
# control for psuedo-high-frequency power from discontinuous (motion-masked) frames
# see https://scipy-cookbook.readthedocs.io/items/ButterworthBandpass.html
from scipy.signal import butter, lfilter
def butter_bandpass(lowcut, highcut, fs, order=2):
# parameters matched to real-data pipeline
nyquist = 0.5 * fs
low = lowcut / nyquist
high = highcut / nyquist
b, a = butter(order, [low, high], btype='band')
return b, a
def butter_bandpass_filter(data, lowcut, highcut, fs, order=5):
b, a = butter_bandpass(lowcut, highcut, fs, order=order)
y = lfilter(b, a, data)
return y
# Desired cutoff frequencies (in Hz).
lowcut = 0.008
highcut = 0.09
# create filter
butter_bandpass(lowcut, highcut, fs)
# filter data
templateDataCort = butter_bandpass_filter(templateDataCort, lowcut, highcut, fs)
# extract frequency distribution
# get frequency distribution of all vertices
freqs, density = signal.welch(templateDataCort,fs,axis=0)
# average frequency distribution
avgFreqs=density.mean(axis=1)
plt.plot(freqs,avgFreqs)
plt.xlabel('Hz')
plt.xlim(-0.01,0.25)
plt.show()
# extract cumulative distribution, bin widths (https://stackoverflow.com/questions/13476807/probability-density-function-from-histogram-in-python-to-fit-another-histrogram)
cum_counts = np.cumsum(avgFreqs)
bin_widths = np.repeat(freqs[1] - freqs[0],cum_counts.shape[0])
x = cum_counts*bin_widths
y = freqs[0:]
# interpolate density of frequency distribution
inverse_density_function = sci.interpolate.interp1d(x, y)
# use interpolated density to sample realistic frequencies for simulated waves
randFreqs = np.zeros(100)
for i in range(len( randFreqs )):
u = random.uniform( x[0], x[-1] )
randFreqs[i] = inverse_density_function( u )
# plot randomly sampled values
plt.hist(randFreqs,100)
plt.xlabel('Randomly Sampled Values from Distr: Hz')
plt.xlim(-0.01,0.25)
plt.show()
vox offset (=1.70334e+06) not divisible by 16, not SPM compatible; leaving at current value vox offset (=1.70334e+06) not divisible by 16, not SPM compatible; leaving at current value C:\Users\adam\anaconda3\lib\site-packages\nibabel\nifti1.py:590: UserWarning: Extension size is not a multiple of 16 bytes; Assuming size is correct and hoping for the best warnings.warn( vox offset (=1.70334e+06) not divisible by 16, not SPM compatible; leaving at current value
